capture log close
log using "${path_log}/02_read_bop_hh.log", replace

*** Compute consumption shares
clear all

forval w = 22/29 {
	
	append using "${path_orig}\BOPHH\bophh_suf_202203_v01_wave`w'.dta" 
	
} // w

recode spentlastmon_? (-9998 -9997 -6666 = .)  

*** Expand panel
xtset id wave
tsfill, full
gsort id wave

*** Spending refers to spending in the past month
replace wave = wave - 1

g spending_total = spentlastmon_a + spentlastmon_b + spentlastmon_c + spentlastmon_d + spentlastmon_e + spentlastmon_f + spentlastmon_g + spentlastmon_h 

*** Calculate average share over sample period
qui foreach v of varlist spentlastmon_a spentlastmon_b spentlastmon_c spentlastmon_d spentlastmon_e spentlastmon_f spentlastmon_g spentlastmon_h {
	
	g share_`v' = `v' / spending_total
 
} // v

collapse (mean) share_*, by(id)

ren (share_spentlastmon_a share_spentlastmon_b share_spentlastmon_c share_spentlastmon_d share_spentlastmon_e share_spentlastmon_f share_spentlastmon_g share_spentlastmon_h) (share_durables share_non_durables share_clothes_shoes share_leisure share_mobility share_services share_vacation share_housing)

tempfile shares 

save `shares'

*** Load waves 21 to 28
use "${path_orig}\BOPHH\bophh_suf_202203_v01_wave21.dta", clear

*** Merge with county IDs
merge 1:1 id wave using "${path_orig}\BOPHH\bophh_mapping_kkz_wave21_v02.dta"
keep if _merge != 2
drop _merge

qui forval w = 22/28 {
	
	append using "${path_orig}\BOPHH\bophh_suf_202203_v01_wave`w'.dta" 
	
	merge 1:1 id wave using "${path_orig}\BOPHH\bophh_mapping_kkz_wave`w'_v01.dta", update
	keep if _merge != 2
	drop _merge
	
} // w

*** Recode missings
recode inflexppoint inflexppoint_long eduschool employ hhsize familystatus hhinc expmacroquali_? netwealth_detail_? expint_sav incexp_? interr homeown devinfpoint infdef (-9999 -9998 -9997 -9996 -6666 = .)  

*** Assign first county households are observed in 
bysort id (wave): g kkz_aux = kkz if _n == 1
gegen county_id = max(kkz_aux), by(id)

*** Expand panel
xtset id wave
tsfill, full
gsort id wave

*** Merge with local labor market IDs
merge m:1 county_id using "${path_data}\county_llm"
keep if _merge != 2
drop _merge  

*** Assign 5-year and 10-year inflation expectations
g d_5y = (drandom1 == 1 & inlist(wave,21,22,25,26,27,28)) | (drandom2 == 1 & inlist(wave,20,23)) | (inlist(drandom1,1,3) & inlist(wave,24))   
g d_10y = (drandom1 == 2 & inlist(wave,21,22,25,26,27,28)) | (drandom2 == 2 & inlist(wave,20,23)) | (inlist(drandom1,2,4) & inlist(wave,24)) 

g inflexppoint5y = inflexppoint_long if d_5y == 1
g inflexppoint10y = inflexppoint_long if d_10y == 1

replace inflexppoint_long = . if d_5y == 0 & d_10y == 0 

*** Merge with consumption shares
merge m:1 id using `shares'
keep if _merge != 2
drop _merge 

*** Merge with CPI sub-indices
merge m:1 wave using "${path_data}\cpi_subindices"
keep if _merge != 2
drop _merge 
 
*** Compute personal inflation perception and expectation for consumption categories
qui foreach s in "durables" "non_durables" "clothes_shoes" "leisure" "mobility" "services" "vacation" "housing" {
	
	g inflation_f12_pers_`s' = share_`s' * inflation_f12_`s'
	g inflation_l12_pers_`s' = share_`s' * inflation_l12_`s'
	
} // s

*** Compute weighted sum of consumption categories
gegen inflation_f12_pers = rowtotal(inflation_f12_pers_*) if !mi(inflation_f12_pers_durables,inflation_f12_pers_non_durables,inflation_f12_pers_clothes_shoes,inflation_f12_pers_leisure,inflation_f12_pers_mobility,inflation_f12_pers_services,inflation_f12_pers_vacation,inflation_f12_pers_housing)
gegen inflation_l12_pers = rowtotal(inflation_l12_pers_*) if !mi(inflation_l12_pers_durables,inflation_l12_pers_non_durables,inflation_l12_pers_clothes_shoes,inflation_l12_pers_leisure,inflation_l12_pers_mobility,inflation_l12_pers_services,inflation_l12_pers_vacation,inflation_l12_pers_housing)
replace inflation_f12_pers = inflation_f12_pers * 100
replace inflation_l12_pers = inflation_l12_pers * 100

*** Merge with inflation data
merge m:1 wave using "${path_data}\cpi_hh"
keep if _merge != 2
drop _merge 

*** Merge with professional forecast data
merge m:1 wave using "${path_data}\forecast_hh"
keep if _merge != 2
drop _merge 

*** Account for negative inflation expectations
replace inflexppoint = - inflexppoint if infdef == 2
g inflexppoint_orig = inflexppoint

*** Account for outliers 
winsor2 inflexppoint5y inflexppoint10y inflexppoint_long expint_sav devinfpoint inflexppoint, cuts(2 98) trim replace  
 
qui foreach v of varlist incexp_? {
	
	replace `v' = `v' / 100
	
} // v

xtset id wave 

*** Appendix Figure A.3a
preserve
	
collapse (mean) inflexppoint, by(wave)

twoway (connected inflexppoint wave, lcolor(blue) mcolor(blue)), graphregion(color(white)) bgcolor(white) plotregion(fcolor(white)) xlabel(21 "Sep 2021" 22 "Oct 2021" 23 "Nov 2021" 24 "Dec 2021" 25 "Jan 2022" 26 "Feb 2022" 27 "Mar 2022" 28 "Apr 2022", angle(90)) xtitle("") ytitle("1-year inflation expectation (in pp)") ylabel(4(1)7)
 
graph export "${path_export_graphs}\appendix_figure_a3a.pdf", replace

restore 
  
*** Merge with wave 28 specific questions
merge m:1 id using "${path_orig}\BOPHH\anonymisation_bophh_originPro_2022_003_wave28_v01"
keep if _merge == 3
drop _merge
 
*** Recode variables
recode exp_energy electricity_price_change reason_change month_price_change (-9998 -9997 -6666 = .)
  
*** No price change equals 0
replace electricity_price_change = 0 if month_price_change == 1 

*** Dummy for no price change
g treated = month_price_change == 1 if !mi(month_price_change)

*** Gender
recode gender interr (2=1) (1=0)

*** East Germany
g east = region == 4 if !mi(region)

*** Education
replace eduschool = 1 if inlist(eduschool,7,8) // current student, no degree or other degree
replace eduschool = 3 if inlist(eduschool,4) // Realschulabschluss + Abschluss Polytechnische Oberschule
g d_educ_high = inlist(eduschool,5,6) if !mi(eduschool) // Abitur + Fachhochschulreife

*** Employment dummies
g fulltime = employ == 1 if !mi(employ)
g parttime = employ == 2 if !mi(employ)
g unemployed = employ == 5 if !mi(employ)
g retired = inlist(employ,7,8) if !mi(employ) // retired or early retirement
g work_other = inlist(employ,3,4,6,9,10,11) if !mi(employ) // Other: internship, parental leave, sick leave, voluntary service, stay-at-home mom/dad, other

*** Homeownership 
g d_homeowner = inlist(homeown,3,4) if !mi(homeown)
ereplace d_homeowner = min(d_homeowner), by(id) // Only households that were homeowners during the whole sample period

*** Marital status
g d_married = inlist(familystatus,1,2) if !mi(familystatus)
ereplace d_married = min(d_married), by(id) // Only households that were married during the whole sample period
 
*** Liquid wealth
g bank_deposits = 0 if netwealth_detail_a == 1
replace bank_deposits = 1250 if netwealth_detail_a == 2
replace bank_deposits = 3750 if netwealth_detail_a == 3
replace bank_deposits = 7500 if netwealth_detail_a == 4
replace bank_deposits = 17500 if netwealth_detail_a == 5
replace bank_deposits = 37500 if netwealth_detail_a == 6
replace bank_deposits = 75000 if netwealth_detail_a == 7
replace bank_deposits = 175000 if netwealth_detail_a == 8 	
replace bank_deposits = 375000 if netwealth_detail_a == 9
replace bank_deposits = 500000 if netwealth_detail_a == 10

g securities = 0 if netwealth_detail_c == 1
replace securities = 1250 if netwealth_detail_c == 2
replace securities = 3750 if netwealth_detail_c == 3
replace securities = 7500 if netwealth_detail_c == 4
replace securities = 17500 if netwealth_detail_c == 5
replace securities = 37500 if netwealth_detail_c == 6
replace securities = 75000 if netwealth_detail_c == 7
replace securities = 175000 if netwealth_detail_c == 8 	
replace securities = 375000 if netwealth_detail_c == 9
replace securities = 500000 if netwealth_detail_c == 10
 
g wealth_liquid = bank_deposits + securities  
gegen wealth_liquid_mean = mean(wealth_liquid), by(id) 
qui sum wealth_liquid_mean, d
g d_low_liquidwealth = wealth_liquid_mean <= `r(p50)' if !mi(wealth_liquid_mean) // median split

*** Education
g d_educ_no = eduschool == 1 if !mi(eduschool)
g d_educ_secondary = eduschool == 2 if !mi(eduschool)
g d_educ_intsecondary = eduschool == 3 if !mi(eduschool)
g d_educ_technical = eduschool == 5 if !mi(eduschool)
g d_educ_highschool = eduschool == 6 if !mi(eduschool)

*** Generate treatment dummies
g price_change = month_price_change == 2 & wave == 21 if !mi(month_price_change)
replace price_change = 1 if month_price_change == 3 & wave == 22
replace price_change = 1 if month_price_change == 4 & wave == 23
replace price_change = 1 if month_price_change == 5 & wave == 24
replace price_change = 1 if month_price_change == 6 & wave == 25
replace price_change = 1 if month_price_change == 7 & wave == 26
replace price_change = 1 if month_price_change == 8 & wave == 27
replace price_change = 1 if month_price_change == 9 & wave == 28
bysort id (wave): g tau = sum(price_change)

*** Income dummies
g d_inc_low = hhinc <= 8 if !mi(hhinc)
g d_inc_b1_2 = inlist(hhinc,3,4) if !mi(hhinc)
g d_inc_b2_3 = inlist(hhinc,5,6) if !mi(hhinc)
g d_inc_b3_4 = inlist(hhinc,7,8) if !mi(hhinc)
g d_inc_b4_5 = inlist(hhinc,9) if !mi(hhinc)
g d_inc_b5_6 = inlist(hhinc,10) if !mi(hhinc)
g d_inc_b6_10 = inlist(hhinc,11,12,13) if !mi(hhinc)

g hhincome = 250 if hhinc == 1
replace hhincome = 750 if hhinc == 2
replace hhincome = 1250 if hhinc == 3
replace hhincome = 1750 if hhinc == 4
replace hhincome = 2250 if hhinc == 5
replace hhincome = 2750 if hhinc == 6
replace hhincome = 3250 if hhinc == 7
replace hhincome = 3750 if hhinc == 8
replace hhincome = 4500 if hhinc == 9
replace hhincome = 5500 if hhinc == 10
replace hhincome = 7000 if hhinc == 11
replace hhincome = 9000 if hhinc == 12
replace hhincome = 10000 if hhinc == 13

*** Prediction errors 
g abs_predict_error = abs(inflexppoint - cpi_f12m)
g sq_predict_error = (inflexppoint - cpi_f12m)^2 
g abs_difference_spf = abs(inflexppoint - consensus_forecast_int)
g sq_difference_spf = (inflexppoint - consensus_forecast_int)^2 
g abs_predict_error_pers = abs(inflexppoint - inflation_f12_pers)
g sq_predict_error_pers = (inflexppoint - inflation_f12_pers)^2

*** Perception errors
g abs_backcast_error = abs(devinfpoint - cpi_l12m)
replace abs_backcast_error = . if tau == 1 // Only use pre-treatment data

gegen mean_abs_backcast_error = mean(abs_backcast_error), by(id) // Calculate pre-treatment mean
qui sum mean_abs_backcast_error, det
g d_informed = mean_abs_backcast_error <= `r(p50)' if !mi(mean_abs_backcast_error) // median split 

*** Perception errors net of time trends & demographics
qui reghdfe abs_backcast_error, absorb(wave) resid
predict abs_backcast_error_w, res

qui reghdfe abs_backcast_error gender age east hhsize d_homeowner d_married d_educ_secondary d_educ_intsecondary d_educ_technical d_educ_highschool fulltime parttime unemployed retired d_inc_b*, absorb(wave) resid
predict abs_backcast_error_d, res

gegen mean_abs_backcast_error_w = mean(abs_backcast_error_w), by(id) // Calculate pre-treatment mean
qui sum mean_abs_backcast_error_w, det
g d_informed_w = mean_abs_backcast_error_w <= `r(p50)' if !mi(mean_abs_backcast_error_w) // median split

gegen mean_abs_backcast_error_d = mean(abs_backcast_error_d), by(id) // Calculate pre-treatment mean
qui sum mean_abs_backcast_error_d, det
g d_informed_d = mean_abs_backcast_error_d <= `r(p50)' if !mi(mean_abs_backcast_error_d) // median split

*** Perception errors compared to personal inflation
g abs_backcast_error_prs = abs(devinfpoint - inflation_l12_pers)  
replace abs_backcast_error_prs = . if tau == 1 // Only use pre-treatment data

gegen mean_abs_backcast_error_prs = mean(abs_backcast_error_prs), by(id)  // Calculate pre-treatment mean
qui sum mean_abs_backcast_error_prs, det
g d_informed_prs = mean_abs_backcast_error_prs <= `r(p50)' if !mi(mean_abs_backcast_error_prs) // median split

*** Dummy for price decreases
g d_price_decrease = inlist(electricity_price_change,1,2) if !mi(electricity_price_change) 

*** Assign treatment wave
g treatment_time = 21 if month_price_change == 2
replace treatment_time = 22 if month_price_change == 3
replace treatment_time = 23 if month_price_change == 4
replace treatment_time = 24 if month_price_change == 5
replace treatment_time = 25 if month_price_change == 6
replace treatment_time = 26 if month_price_change == 7
replace treatment_time = 27 if month_price_change == 8
replace treatment_time = 28 if month_price_change == 9

*** Assign price increase/decrease 
g price_increase = 0 if electricity_price_change == 0
replace price_increase = -0.10 if electricity_price_change == 1
replace price_increase = -0.05 if electricity_price_change == 2
replace price_increase = 0.05 if electricity_price_change == 3
replace price_increase = 0.15 if electricity_price_change == 4
replace price_increase = 0.25 if electricity_price_change == 5
replace price_increase = 0.35 if electricity_price_change == 6
replace price_increase = 0.45 if electricity_price_change == 7
replace price_increase = 0.50 if electricity_price_change == 8

*** Alternative assumption about endpoints: +- 5
g price_increase5 = price_increase
replace price_increase5 = -0.15 if electricity_price_change == 1
replace price_increase5 =  0.55 if electricity_price_change == 8

*** Alternative assumption about endpoints: +- 10
g price_increase10 = price_increase
replace price_increase10 = -0.20 if electricity_price_change == 1
replace price_increase10 =  0.60 if electricity_price_change == 8

replace price_increase = price_increase * 100 
replace price_increase5 = price_increase5 * 100 
replace price_increase10 = price_increase10 * 100 

*** Low vs high increases
qui sum price_increase if price_increase > 0, det
g d_low_increase = inlist(electricity_price_change,0) | price_increase <= `r(p50)' if !mi(price_increase) // median split
g d_high_increase = inlist(electricity_price_change,0) | price_increase > `r(p50)' if !mi(price_increase) // median split
	
*** Scaled Treatment Indicator
g tau_price_increase = tau * price_increase   
 
*** Dummy for high energy cost share
g share_exp_energy = exp_energy / hhincome 
winsor2 share_exp_energy, cuts(0 98) replace // Account for outliers

qui sum share_exp_energy, d
g d_share_exp_energy_high = share_exp_energy > `r(p50)' if !mi(share_exp_energy) // median split

*** Adjust scaling
replace share_exp_energy = share_exp_energy * 100
replace hhincome = hhincome / 1000
replace wealth_liquid_mean = wealth_liquid_mean / 1000
replace exp_energy = exp_energy / 100
 
*** Compute fixed effects 
g d_inc_high = 1 - d_inc_low
gegen d_inc_low_wave = group(d_inc_low wave) 
g d_educ_low = 1 - d_educ_high
gegen d_educ_high_wave = group(d_educ_high wave)
g d_high_liquidwealth = 1 - d_low_liquidwealth
gegen d_low_liquidwealth_wave = group(d_low_liquidwealth wave)
g d_uninformed = 1 - d_informed
gegen d_informed_wave = group(d_informed wave) 
g d_uninformed_prs = 1 - d_informed_prs
gegen d_informed_prs_wave = group(d_informed_prs wave) 
g d_uninformed_w = 1 - d_informed_w
gegen d_informed_w_wave = group(d_informed_w wave) 
g d_uninformed_d = 1 - d_informed_d
gegen d_informed_d_wave = group(d_informed_d wave) 
g d_renter = 1 - d_homeowner
gegen d_homeowner_wave = group(d_homeowner wave)
gegen d_share_exp_energy_high_wave = group(d_share_exp_energy_high wave)
 
xtset id wave 

*** Drop HHs observed more than three months after treatment
g price_change_L3 = L3.price_change
replace price_change_L3 = 0 if price_change_L3 == .
bysort id: g price_change_L3_tmp = sum(price_change_L3)
replace price_change_L3 = price_change_L3_tmp  
drop if price_change_L3 == 1  
	 
gsort id wave

*** Drop missing observations 
drop if mi(gender) 

*** Keep only panel HHs
bysort id: g num_waves = _N
keep if num_waves != 1
 
*** Compute implied mean and standard deviation of income expectations following Coibion et al (2024) "The Effect of Macroeconomic Uncertainty on Household Spending", American Economic Review vol. 114, no. 3, March 2024 (pp. 645–77)
preserve

keep id wave incexp_?

gegen check = rowtotal(incexp_*)
replace check = round(check,0.01)
keep if check == 1 

reshape long incexp_, i(id wave) j(bin) string

ren incexp_ incexp

g lb = .
g ub = .

replace lb = -2500 if bin == "a"
replace ub = -2000 if bin == "a"
replace lb = -2000 if bin == "b"
replace ub = -1500 if bin == "b"
replace lb = -1500 if bin == "c"
replace ub = -1000 if bin == "c"
replace lb = -1000 if bin == "d"
replace ub = -500 if bin == "d"
replace lb = -500 if bin == "e"
replace ub = -250 if bin == "e"
replace lb = -250 if bin == "f"
replace ub = 0 if bin == "f"
replace lb = 0 if bin == "g"
replace ub = 250 if bin == "g"
replace lb = 250 if bin == "h"
replace ub = 500 if bin == "h"
replace lb = 500 if bin == "i"
replace ub = 1000 if bin == "i"
replace lb = 1000 if bin == "j"
replace ub = 1500 if bin == "j" 
replace lb = 1500 if bin == "k"
replace ub = 2000 if bin == "k" 
replace lb = 2000 if bin == "l"
replace ub = 2500 if bin == "l" 	 

gsort id wave lb
gegen num_bin = group(lb) 

g incexp_mean = .
g incexp_var = .

gegen id_wave = group(id wave)
tsset id_wave num_bin
 
g bin_pos = incexp != 0 if !mi(incexp)
bysort id_wave (num_bin): gegen num_bin_pos = sum(bin_pos)
bysort id_wave (num_bin): g num_bin_cumul = sum(bin_pos)

*** Cumulative density distribution
g cdf0 = incexp if num_bin == 1
replace cdf0 = L.cdf0 + incexp if num_bin > 1 

*** One bin
replace incexp_mean = lb + (ub - lb)/2 if incexp == 1
replace incexp_var = (ub - lb)^2/24 if incexp == 1

*** Two bins 
g temp_incexp = incexp if num_bin_pos == 2
replace temp_incexp = . if temp_incexp == 0

bysort id_wave: gegen temp_incexp_min = min(temp_incexp)
replace temp_incexp_min = . if mi(temp_incexp)

*** Equal probabilities: take an union of invervals and then the midpoint as the implied mean
bysort id_wave: gegen temp_incexp_min1 = min(lb) if temp_incexp_min == 0.5
bysort id_wave: gegen temp_incexp_max1 = max(ub) if temp_incexp_min == 0.5
 
replace incexp_mean = temp_incexp_min1 + 0.5*(temp_incexp_max1 - temp_incexp_min1) if temp_incexp_min == 0.5 & num_bin_pos == 2
replace incexp_var = (temp_incexp_max1 - temp_incexp_min1)^2/24 if temp_incexp_min == 0.5 & num_bin_pos == 2
	
drop temp_incexp_min1 temp_incexp_max1

*** Unequal probabilities: the "high" bin has a larger probability   
g temp_flag = incexp > 0.5 if num_bin_pos == 2 & num_bin_cumul == 2 & !mi(incexp)
g temp_incexp_mass = 1 - incexp if temp_flag == 1

g m = (lb - ub*sqrt(temp_incexp_mass/2))/(1 - sqrt(temp_incexp_mass/2)) if temp_flag == 1
replace incexp_mean = m + 0.5*(ub - m) if temp_flag == 1
replace incexp_var = (m - ub)^2/24 if temp_flag == 1

capture drop temp_flag
capture drop temp_incexp_mass
capture drop m

*** Unequal probabilities: the "low" bin has a larger probability 
g temp_flag = incexp > 0.5 if num_bin_pos == 2 & num_bin_cumul == 1 & !mi(incexp)
g temp_incexp_mass = 1 - incexp if temp_flag == 1

g m = (ub - lb*sqrt(temp_incexp_mass/2))/(1 - sqrt(temp_incexp_mass/2)) if temp_flag == 1
replace incexp_mean = lb + 0.5*(m - lb) if temp_flag == 1
replace incexp_var = (m - lb)^2/24 if temp_flag == 1
 
capture drop temp_flag
capture drop temp_incexp_mass
capture drop m
capture drop temp_incexp*
  
*** Three or more bins
capture drop temp_incexp*
g temp_incexp = incexp if num_bin_pos >= 3 & !mi(num_bin_pos)
replace temp_incexp = . if temp_incexp == 0

bysort id_wave: gegen temp_incexp_min = min(temp_incexp)
replace temp_incexp_min = . if mi(temp_incexp)

bysort id_wave: gegen temp_incexp_min_lb = min(lb) if !mi(temp_incexp) 
bysort id_wave: gegen temp_incexp_max_ub = max(ub) if !mi(temp_incexp) 

g t = (ub - temp_incexp_min_lb) / (temp_incexp_max_ub - temp_incexp_min_lb)

gegen id_wave_num = group(id_wave) if num_bin_pos >= 3 & !mi(num_bin_pos)

g a1 = .
g b1 = . 

g A = . 
g B = . 

qui sum id_wave_num
local max_obs = `r(max)'

*** Loop over observations to estimate beta distribution for each observation
qui forvalues i=1(1)`r(max)' {
	 
	noi di "obs `i' of `max_obs'"
	nl (cdf0 = ibeta({a},{b},t)) if id_wave_num == `i' & !mi(t), initial(a 1.23 b 2)

	matrix M = e(b)
	replace a1 = M[1,1] if e(sample)
	replace b1 = M[1,2] if e(sample)
	sum temp_incexp_min_lb if e(sample)
	replace A = r(mean) if e(sample)
	sum temp_incexp_max_ub if e(sample)
	replace B = r(mean) if e(sample)

	replace incexp_mean = a1/(a1 + b1)*(B-A) + A if e(sample) // mean of beta distribution
	replace incexp_var = (a1 * b1)/((a1 + b1)^2 * (a1 + b1 + 1))*(B - A)^2 if e(sample) // variance of beta distribution
 
}
 
g incexp_std = sqrt(incexp_var) // get standard deviation
 
collapse (mean) incexp_mean incexp_std, by(id wave)

save "${path_data}\incexp.dta", replace 

restore 

*** Merge with income expectations 
merge 1:1 id wave using "${path_data}\incexp.dta" 
keep if _merge != 2
drop _merge

*** Express in 100€
replace incexp_mean = incexp_mean / 100
replace incexp_std = incexp_std / 100 

qui foreach v of varlist share_durables share_non_durables share_clothes_shoes share_leisure share_mobility share_services share_vacation {
		
	replace `v' = `v' * 100
	
}

*** Dummies for liquid wealth
qui sum wealth_liquid_mean, det
g d_wealth_liquid_p2 = wealth_liquid_mean >= `r(p25)' & wealth_liquid_mean < `r(p50)' if !mi(wealth_liquid_mean)
g d_wealth_liquid_p3 = wealth_liquid_mean >= `r(p50)' & wealth_liquid_mean < `r(p75)' if !mi(wealth_liquid_mean)
g d_wealth_liquid_p4 = wealth_liquid_mean >= `r(p75)' if !mi(wealth_liquid_mean)

*** Generate Huber weights
qui rreg inflexppoint_orig tau if d_price_decrease != 1, genwt(wgt_outlier)

*** Number of observations within each local labor market
qui levelsof llm_id, clean local(llm_list)
g num_llm = .
qui foreach llm of local llm_list {
	
	gdistinct id if llm_id == `llm'
	replace num_llm = `r(ndistinct)' if llm_id == `llm'
	
} // a 

*** Require at least 10 observations in each local labor market 
replace llm_id = . if num_llm < 10
gegen llm_wave = group(llm_id wave)

*** Save data set including HHs that reported a price decrease
preserve 

bysort id (wave): g person_time = _n
xtset id person_time 

keep id wave treatment_time person_time price_increase d_price_decrease tau inflexppoint inflexppoint_orig wgt_outlier reason_change interr

save "${path_data}\bophh_wdec.dta", replace

restore 

*** Drop households reporting price decreases
drop if d_price_decrease == 1

*** Generate bins for personal inflation distribution
g inflation_l12_pers_d = floor(inflation_l12_pers / 0.25)

gegen aux_distinct = tag(id inflation_l12_pers_d) if !mi(inflation_l12_pers_d)
gegen distinct = total(aux_distinct) if !mi(inflation_l12_pers_d), by(inflation_l12_pers_d) 
  
*** Appendix Figure A.4		
hist inflation_l12_pers if !(distinct < 5), start(2) width(0.25) color(blue%50) percent graphregion(color(white)) bgcolor(white) plotregion(fcolor(white)) xtitle(" " "Personal inflation rate (in pp)") xlabel(2(2)12) nodraw
drop aux_distinct distinct

graph display, xsize(20) ysize(11)

graph export "${path_export_graphs}\appendix_figure_a4.pdf", replace

*** Appendix Table A.4
qui estpost sum inflexppoint inflexppoint5y inflexppoint10y expint_sav incexp_mean incexp_std expmacroquali_f expmacroquali_a expmacroquali_g expmacroquali_i expmacroquali_b expmacroquali_x abs_predict_error sq_predict_error abs_predict_error_pers sq_predict_error_pers abs_difference_spf sq_difference_spf hhincome wealth_liquid_mean exp_energy, det   
est store a
 
esttab a using "${path_export_tables}\appendix_table_a4.tex", replace nomtitles collabels(\multicolumn{1}{c}{mean} \multicolumn{1}{c}{median} \multicolumn{1}{c}{sd} \multicolumn{1}{c}{N}) cells("mean(fmt(2)) p50(fmt(2)) sd(fmt(2)) count(fmt(0))") coeflabels(inflexppoint "1-year inflation expectation (in pp)" inflexppoint5y "5-year inflation expectation (in pp)" inflexppoint10y "10-year inflation expectation (in pp)" expint_sav "interest rate expectation (in pp)" incexp_mean "income change expectation (in 100\euro)" incexp_std "income change expectation uncertainty (in 100\euro)" expmacroquali_a "unemployment rate expectation" expmacroquali_f "house price expectation" expmacroquali_g "growth rate expectation" expmacroquali_i "stock market expectation" expmacroquali_b "rent expectation" expmacroquali_x "taxation expectation" hhincome "household net income (in 1000\euro)" exp_energy "monthly electricity expenditures (in 100\euro)" wealth_liquid_mean "household liquid wealth (in 1000\euro)" abs_predict_error "abs. forecast error (in pp)" sq_predict_error "sq. forecast error (in pp)" abs_difference_spf "abs. difference to professional forecasters (in pp)" sq_difference_spf "sq. difference to professional forecasters (in pp)" abs_predict_error_pers "abs. forecast error: personal inflation (in pp)" sq_predict_error_pers "sq. forecast error: personal inflation (in pp)") nonumber f noobs alignment(S) booktabs postfoot("\bottomrule \bottomrule")

est clear

*** Figure 1a
preserve 

collapse (max) month_price_change, by(id)

hist month_price_change, percent xlabel(1 "no price change" 2 "Sep 2021" 3 "Oct 2021" 4 "Nov 2021" 5 "Dec 2021" 6 "Jan 2022" 7 "Feb 2022" 8 "Mar 2022" 9 "Apr 2022", angle(90)) xtitle("") color(blue) discrete gap(50) graphregion(color(white)) bgcolor(white) plotregion(fcolor(white))

graph export "${path_export_graphs}\figure_1a.pdf", replace 

restore

*** Appendix Figure A.5a
hist wave, frequency color(blue) graphregion(color(white)) bgcolor(white) plotregion(fcolor(white)) xlabel(21 "Sep 2021" 22 "Oct 2021" 23 "Nov 2021" 24 "Dec 2021" 25 "Jan 2022" 26 "Feb 2022" 27 "Mar 2022" 28 "Apr 2022", angle(90)) xtitle("") ytitle("number of observations" " ") discrete gap(50) 

graph export "${path_export_graphs}\appendix_figure_a5a.pdf", replace 

*** Normalize by standard devation
qui foreach v of varlist expmacroquali_? { 
	
	sum `v'
	replace `v' = `v' / `r(sd)'
	
}

xtset id wave 

save "${path_data}\bophh_final.dta", replace

log close